Debates—Stochastic subsurface hydrology from theory to practice: why stochastic modeling has not yet permeated into practitioners?

We address modern topics of stochastic hydrogeology from their potential relevance to real modeling efforts at the field scale. While the topics of stochastic hydrogeology and numerical modeling have become routine in hydrogeological studies, nondeterministic models have not yet permeated into practitioners. We point out a number of limitations of stochastic modeling when applied to real applications and comment on the reasons why stochastic models fail to become an attractive alternative for practitioners. We specifically separate issues corresponding to flow, conservative transport, and reactive transport. The different topics addressed are emphasis on process modeling, need for upscaling parameters and governing equations, relevance of properly accounting for detailed geological architecture in hydrogeological modeling, and specific challenges of reactive transport. We end up by concluding that the main responsible for nondeterministic models having not yet permeated in industry can be fully attributed to researchers in stochastic hydrogeology

Mathematical equivalence between time-dependent single-rate and multirate mass transfer models

The often observed tailing of tracer breakthrough curves is caused by a multitude of mass transfer processes taking place over multiple scales. Yet, in some cases, it is convenient to fit a transport model with a single-rate mass transfer coefficient that lumps all the non-Fickian observed behavior. Since mass transfer processes take place at all characteristic times, the single-rate mass transfer coefficient derived from measurements in the laboratory or in the field vary with time w(t). The literature review and tracer experiments compiled by Haggerty et al. (2004) from a number of sites worldwide suggest that the characteristic mass transfer time, which is proportional to w(t)^-1, scales as a power law of the advective and experiment duration. This paper studies the mathematical equivalence between the multirate mass transfer model (MRMT) and a time-dependent single-rate mass transfer model (t-SRMT). In doing this, we provide new insights into the previously observed scale-dependence of mass transfer coefficients. The memory function, g(t), which is the most salient feature of the MRMT model, determines the influence of the past values of concentrations on its present state. We found that the t-SRMT model can also be expressed by means of a memory function \phi(t,\tau). In this case, though the memory function is nonstationary, meaning that in general it cannot be written as \phi(t-\tau). Nevertheless, the full behavior of the concentrations using a single time-dependent rate w(t) is approximately analogous to that of the MRMT model provided that the equality w(t) = -dlng(t)/dt holds and the field capacity is properly chosen. This relationship suggests that when the memory function is a power law, g(t) \approx t^{1-k}, the equivalent mass transfer coefficient scales as w(t) \approx t^-1, nicely fitting without calibration the estimated mass transfer coefficients compiled by Haggerty et al. (2004)

Optimal reconstruction of concentrations, gradients and reaction rates from particle distributions

Random walk particle tracking methodologies to simulate solute transport of conservative species constitute an attractive alternative for their computational efficiency and absence of numerical dispersion. Yet, problems stemming from the reconstruction of concentrations from particle distributions have typically prevented its use in reactive transport problems. The numerical problem mainly arises from the need to first reconstruct the concentrations of species/components from a discrete number of particles, which is an error prone process, and then computing a spatial functional of the concentrations and/or its derivatives (either spatial or temporal). Errors are then propagated, so that common strategies to reconstruct this functional require an unfeasible amount of particles when dealing with nonlinear reactive transport problems. In this context, this article presents a methodology to directly reconstruct this functional based on kernel density estimators. The methodology mitigates the error propagation in the evaluation of the functional by avoiding the prior estimation of the actual concentrations of species. The multivariate kernel associated with the corresponding functional depends on the size of the support volume, which defines the area over which a given particle can influence the functional. The shape of the kernel functions and the size of the support volume determines the degree of smoothing, which is optimized to obtain the best unbiased predictor of the functional using an iterative plug-in support volume selector. We applied the methodology to directly reconstruct the reaction rates of a precipitation/dissolution problem involving the mixing of two different waters carrying two aqueous species in chemical equilibrium and moving through a randomly heterogeneous porous mediu

Convergent-flow tracer tests in heterogeneous media: combined experimental–numerical analysis for determination of equivalent transport parameters

In modeling transport within naturally heterogeneous aquifers, it is usually assumed that the transport equations valid at local scales can also be applied at larger scales. At larger scales, the heterogeneous domain is represented by an equivalent homogeneous medium. Convergent-flow tracer tests constitute one of the most frequently used field tests to estimate effective input parameters of equivalent homogeneous aquifers. Traditionally, statistical approaches applied to groundwater flow and solute transport have provided tools to estimate these equivalent parameters. These approaches are based on a number of simplifications including the assumption that the point transmissivity values follow a multilog-normal random function. Several investigators have found that this assumption may not be valid in many field cases. In order to study the applicability of the equivalent homogeneous formulation in a nontraditional stochastic field, a number of experimental and numerical studies were conducted. The results are used to determine the apparent values of porosity and dispersivity that would be obtained if convergent-flow tracer tests were conducted in a deterministically generated heterogeneous transmissivity field displaying anisotropy in the correlation structure. It is shown that in this particular heterogeneous media, apparent porosity strongly depends on connectivity rather than on transmissivity. This dependence on connectivity questions the theoretical results obtained in continuum equivalent fields to estimate effective porosity

Point-to-point connectivity, an abstract concept or a key issue for risk assessment studies?

Connectivity of high/low-permeability areas has been recognized to significantly impact groundwater flow and solute transport. The task of defining a rigorous quantitative measure of connectivity for continuous variables has failed so far, and thus there exist a suite of connectivity indicators which are dependent on the specific hydrodynamic processes and the interpretation method. Amongst the many existing indicators, we concentrate on those characterizing connectivity between the points involved in a hydraulic or tracer test. The flow connectivity indicator used here is based on the time elapsed for hydraulic response in a pumping test (e.g., the storage coefficient estimated by the Cooper–Jacob method, Sest). Regarding transport, we select the estimated porosity from the breakthrough curve (ϕest). According to Knudby and Carrera [Knudby C, Carrera J. On the relationship between indicators of geostatistical, flow and transport connectivity. Adv Water Resour 2005;28(4):405–21] these two indicators measure connectivity differently, and are poorly correlated. Here, we use perturbation theory to analytically investigate the intrinsic relationship between Sest and ϕest. We find that ϕest can be expressed as a weighted line integral along the particle trajectory involving two parameters: the transmissivity point values, T, and the estimated values of Sest along the particle path. The weighting function is linear with the distance from the pumping well, thus the influence of the weighting function is maximum at the injection area, whereas the hydraulic information close to the pumping well becomes redundant (null weight). The relative importance of these two factors is explored using numerical simulations in a given synthetic aquifer and tested against intermediate-scale laboratory tracer experiments. We conclude that the degree of connectivity between two points of an aquifer (point-to-point connectivity) is a key issue for risk assessment studies aimed at predicting the travel time of a potential contaminant

Interpretation of column experiments of transport of solutes undergoing an irreversible bimolecular reaction using a continuum approximation

We provide a quantitative interpretation of the column experiment reported by Gramling et al. (2002). The experiment involves advection‐dominated transport in porous media of three dissolved species, i.e., two reactants undergoing a fast irreversible reaction and the resulting product. The authors found that their observations could not be properly fitted with a model based on an advection‐dispersion‐reaction equation (ADRE) assuming the reaction was instantaneous, the actual measured total reaction product being lower than predictions for all times. The data have been recently well reproduced by Edery et al. (2009, 2010) by means of a particle tracking approach in a continuous time random walk framework. These and other authors have questioned the use of partial differential equation (PDE)–based approaches to quantify reactive transport because of the difficulty in capturing local‐scale mixing and reaction. We take precisely this approach and interpret the experiments mentioned by means of a continuum‐scale model based on the ADRE. Our approach differs from previous modeling attempts in that we imbue effects of incomplete mixing at the pore scale in a time‐dependent kinetic reaction term and show that this model allows quantitative interpretation of the experiments in terms of both reaction product profiles and time‐dependent global production rate. The time dependence of the kinetic term presented accounts for the progressive effects of incomplete mixing due to pore‐scale rate‐limited mass transfer, and follows a power law, which is consistent with the compilation of existing experiments reported by Haggerty et al. (2004). Our interpretation can form the basis for further research to assess the potential use of PDE approaches for the interpretation of reactive transport problems in moderately heterogeneous medi

Conditional stochastic mapping of transport connectivity

We present a method for the stochastic simulation of point‐to‐point transport connectivity honoring data from three types of information: (1) travel time estimates obtained from field tracer tests; (2) estimates of flow connectivity indicators obtained from the relatively fast or slow flow response that is observed at a point location given the flow impulse at another location, and (3) measurements of transmissivity at a local scale. The method thus efficiently integrates data obtained from different hydraulic tests, each sampling different areas within the aquifer. To achieve this, we first extend the concept of point‐to‐point flow connectivity and transport connectivity, mathematically formulated by Trinchero et al. (2008) for pumping conditions, to support a more general flow configuration. Interestingly, point‐to‐point flow connectivity can be generally seen as a weighted integral of transmissivity over the entire domain, the weighting function being proportional to the sensitivity of heads with respect to the natural log of transmissivity per unit of aquifer volume. On the contrary, point‐to‐point transport connectivity is a weighted integral along the particle path of the solute mass that involves two variables: transmissivity and flow connectivity. Each variable has its own distinct weighting function. The weighting function of transmissivity is inversely proportional to both the homogeneous travel time and the point velocity sampled along the travel path. On this basis, we show how to generate conditional point‐to‐point transport connectivity maps. The method avoids the inference of cross‐covariance functions between variables measured over different scales and sampled areas (which cannot be otherwise estimated with a few data measurements) by expressing them as a function of the local transmissivity covariance function. An example of the method is provided to evaluate the worth of including tracer data to delineate capture zones of abstraction wells originally defined from local transmissivity measurements. Monte Carlo simulations reveal that the impact of including tracer data is a maximum when the travel time data are obtained at a location different than that of transmissivity measurements. The reason is that weighting functions give larger weights to the injection location, so introducing tracer test data at points where transmissivity is already known is somewhat redundan

Stochastic estimation of hydraulic transmissivity fields using flow connectivity indicator data

Most methods for hydraulic test interpretation rely on a number of simplified assumptions regarding the homogeneity and isotropy of the underlying porous media. This way, the actual heterogeneity of any natural parameter, such as transmissivity ( math formula), is transferred to the corresponding estimates in a way heavily dependent on the interpretation method used. An example is a long-term pumping test interpreted by means of the Cooper-Jacob method, which implicitly assumes a homogeneous isotropic confined aquifer. The estimates obtained from this method are not local values, but still have a clear physical meaning; the estimated math formula represents a regional-scale effective value, while the log-ratio of the normalized estimated storage coefficient, indicated by math formula, is an indicator of flow connectivity, representative of the scale given by the distance between the pumping and the observation wells. In this work we propose a methodology to use math formula, together with sampled local measurements of transmissivity at selected points, to map the expected value of local math formula values using a technique based on cokriging. Since the interpolation involves two variables measured at different support scales, a critical point is the estimation of the covariance and crosscovariance matrices. The method is applied to a synthetic field displaying statistical anisotropy, showing that the inclusion of connectivity indicators in the estimation method provide maps that effectively display preferential flow pathways, with direct consequences in solute transport

Reaction rates and effective parameters in stratified aquifers

Chemical species are advected by water and undergo mixing processes due to effects of local diffusion and/or dispersion. In turn, mixing causes reactions to take place so that the system can locally equilibrate. In general, a multicomponent reactive transport problem is described through a system of coupled non-linear partial differential equations. Under instantaneous chemical equilibrium, a complex geochemical problem can be highly simplified by fully defining the system in terms of conservative quantities, termed master species or components, and the space–time distribution of reaction rates. We investigate the parameters controlling reaction rates in a heterogeneous aquifer at short distances from the source. Hydraulic conductivity at this scale is modeled as a random process with highly anisotropic correlation structure. In the limit for very large horizontal integral scales, the medium can be considered as stratified. Upon modeling transport by means of an ADE (Advection Dispersion Equation), we derive closed-form analytical solutions for statistical moments of reaction rates for the particular case of negligible transverse dispersion. This allows obtaining an expression for an effective hydraulic conductivity, , as a representative parameter describing the mean behavior of the reactive system. The resulting is significantly smaller than the effective conductivity representative of the flow problem. Finally, we analyze numerically the effect of accounting for transverse local dispersion. We show that transverse dispersion causes no variation in the distribution of (ensemble) moments of local reaction rates at very short travel times, while it becomes the dominant effect for intermediate to large travel times

A Bayesian approach to integrate temporal data into probabilistic risk analysis of monitored NAPL remediation

Upon their release into the subsurface, non-aqueous phase liquids (NAPLs) dissolve slowly in groundwater and/or volatilize in the vadose zone threatening the environment and public health over extended periods of time. The failure of a treatment technology at any given site is often due to the unnoticed presence of dissolved NAPL trapped in low permeability areas and/or the remaining presence of substantial amounts of pure phase NAPL after remediation efforts. The design of remediation strategies and the determination of remediation endpoints are traditionally carried out without quantifying risks associated with the failure of such efforts. We conduct a probabilistic risk analysis (PRA) to estimate the likelihood of failure of an on-site NAPL treatment technology. The PRA integrates all aspects of the problem (causes, pathways, and receptors) without resorting to extensive modeling. It accounts for a combination of multiple mechanisms of failure of a monitoring system, such as bypassing, insufficient sampling frequency and malfunctioning of the observation wells. We use a Bayesian framework to update the likelihood of failure of the treatment technology with observed measurements of concentrations at nearby monitoring well